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Within the framework of probability distributions on projective Hilbert space a scheme for the 
calculation of multitime correlation functions is developed. The starting point is the Markovian 
stochastic wave function description of an open quantum system coupled to an environment consist- 
ing of an ensemble of harmonic oscillators in arbitrary pure or mixed states. It is shown that matrix 
elements of reduced Heisenberg picture operators and general time-ordered correlation functions 
can be expressed by time-symmetric expectation values of extended operators in a doubled Hilbert 
space. This representation allows the construction of a stochastic process in the doubled Hilbert 
space which enables the determination of arbitrary matrix elements and correlation functions. The 
numerical efficiency of the resulting stochastic simulation algorithm is investigated and compared 
with an alternative Monte Carlo wave function method proposed first by Dalibard et al. [Phys. 
Rev. Lett. 68, 580 (1992)]. By means of a standard example the suggested algorithm is shown to 
be more efficient numerically and to converge faster. Finally, some specific examples from quantum 
optics are presented in order to illustrate the proposed method, such as the coupling of a system to 
a vacuum, a squeezed vacuum within a finite solid angle, and a thermal mixture of coherent states. 
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I. INTRODUCTION 

In recent years several stochastic wave function meth- 
ods have been developed and used to describe the dynam- 
ics of open quantum systems jl]-|| . All these approaches 
are based on the following idea: instead of solving the 
quantum master equation to obtain the time evolution of 
the reduced density matrix, an ensemble of pure states 
is propagated using a stochastic time evolution. This 
method provides two major advantages: the individual 
sample paths of the different realizations can - in some 
situations - be interpreted as the time evolution of an in- 
dividual continuously monitored quantum system |^ ^] 
and the scaling of the numerical performance with the 
system size is better for this approach . 

A particular interesting situation considered in quan- 
tum optics is the coupling of a system (e. g., an atom or a 
cavity) to the continuum of modes of the electromagnetic 
field. Since a lot of theoretical and experimental effort is 
used to prepare the environment in certain well defined 
states which are not restricted to thermal mixtures (see, 
e. g., 1 15 17)), it is interesting to investigate the time evo- 
lution of an open quantum system which is coupled to an 
environment in an arbitrary state. One way to obtain the 
stochastic time evolution of the reduced system takes as 
its starting point the quantum master equation for the 
reduced density matrix. A stochastic wave function is 
constructed on a phenomenological basis in the following 
way: start with an ensemble of pure states representing 
the initial density matrix p(to) = po and propagate each 
member of the ensemble using a stochastic time evolu- 
tion which guarantees that the covariance matrix of the 
stochastic process is the reduced density matrix. This 



procedure ensures that expectation values of system ob- 
servables are calculated correctly. 

In applications of open quantum systems one is also 
interested in multitime correlation functions such as 
(A(t + r)B(t)), where A(t + r) and B(t) are Heisenberg 
operators, and the angular brackets denote the ensemble 
average. In the density matrix approach these quantities 
are usually determined by making use of the quantum 
regression theorem |l^,[l!| which is based on the identity 

(A(t + r)B(t)) 



AV(t + T,t){BV(t,t ){po}}\, (1) 



Tr 

-Lis 



where Tr sys denotes the trace over the system's degrees 
of freedom and V(t,to) is the time evolution super op- 
erator of the corresponding quantum master equation, 
and A and B are Schrodinger operators. Eq. (|l|) al- 
lows the following interpretation: for the calculation of 
(A(t + r)B(t)) start with the density matrix p(to) = po 
and propagate to p(t) using the quantum master equa- 
tion. Then propagate the "density matrix" Bp(t) up to 
the time t + t and calculate the quantum mechanical ex- 
pectation value of A with respect to this density matrix. 
However, since Bp(t) is, in general, not a positive ma- 
trix (and hence not a true density matrix), it can not 
be expressed as an ensemble of pure states and hence a 
direct generalization of this computational scheme to the 
stochastic wave function approach is not possible. In- 
stead, some alternative computational schemes are pro- 
posed in the literature such as expressing time-ordered 
multitime correlation functions as the sum of symmet- 
ric time-ordered correlation functions Q| which can be 
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simulated directly, or using the Heisenberg picture of the 
quantum state diffusion model f22|| . 

In this article we will present a different approach to 
the derivation of the stochastic time evolution which is 
formulated within the framework of probability distri- 
butions on projective Hilbert space |Q|§- The starting 
point for this approach is not the quantum Master equa- 
tion for the reduced density matrix, but a microscopic 
model for the dynamics of the total system. This model 
is used to calculate the unitary time evolution of the to- 
tal system in second order perturbation theory, which 
leads to a Liouville equation for the corresponding prob- 
ability distribution. The equation of motion of the re- 
duced system's probability distribution is obtained by 
invoking a specific reduction formula, which relies upon 
the above stated condition that expectations of system 
operators are calculated accurately, and performing the 
Markov approximation. This scheme allows the deriva- 
tion of the stochastic time evolution of an open quantum 
system coupled to an environment which is in an almost 
arbitrary state. To be more precise, the only restrictions 
we impose on the state of the environment concern the 
relative time scales of the system's and the environment's 
time evolution. Note, that this restriction is not funda- 
mental but necessary in order to end up with a Markovian 
dynamics of the reduced system. 

To construct a stochastic process that allows also the 
determination of multitime correlation functions we shall 
follow in this article the following strategy: multitime 
correlation functions are expressed in terms of matrix 
elements of reduced Heisenberg picture operators which 
can be written as expectation values in a doubled Hilbert 
space. This then enables the formulation of a stochas- 
tic process in the doubled Hilbert space which fulfills the 
condition, that matrix elements of reduced Heisenberg 
picture operators and hence multitime correlation func- 
tions are calculated correctly. 

This article is organized as follows: In Sect. || we will 
briefly review the statistical description of quantum sys- 
tems in terms of probability distributions on the under- 



lying Hilbert space. Sec. [II describes the derivation of 



the equation of motion of the reduced system's probabil- 
ity density functional. We find that the time evolution of 
this probability density is governed by a Liouville-Master 
equation. The resulting stochastic process is thus a piece- 
wise deterministic Markov process. In Sec. |^ we define 
the matrix elements of the reduced Heisenberg picture 
operators. We show that these matrix elements can be 
determined by simulating a stochastic process in a dou- 
bled Hilbert space which is again a piecewise determinis- 
tic Markov process. The process defined in the doubled 
Hilbert space is used to calculate arbitrary measurable 
multitime correlation functions, which is shown in Sec. 
In order to illustrate the general theory presented in this 
article we explicitly discuss some examples of quantum 
optical systems in Sec. VI. 



II. PROBABILITY DISTRIBUTIONS ON 
HILBERT SPACE 



The basis for a description of open quantum systems 
in terms of a stochastic wave function is the definition of 
a probability measure on the underlying Hilbert space 
TL. In this article we will make use of the probabil- 
ity density functional P[4>, t) which is defined such that 
P[ip,t]Dif;Dij}* is the probability of finding the state of 
the system in the volume element DipDip* around ip at 
time t The volume element of the Hilbert space is 

defined as 



(2) 



where x labels a complete set of quantum numbers. Thus, 
the normalization of P[4>, t] reads 



DipD^*P[il),t} = 1, 



(3) 



where the integral extends over the Hilbert space TL. In 
order to be consistent with the general principles of quan- 
tum mechanics we have to impose two further conditions 
on P[ijj,t]: i) since two state vectors which only differ 
by a phase factor describe the same physical state, we 
require that P[4>, t] does not depend on the phase of the 
state vector and ii) we are only interested in normalized 
state vectors, so that the support of P[ijj,t] is the unit 
sphere in the Hilbert space TL. The quantum mechanical 
expectation value of an arbitrary linear operator A with 
respect to P[ip,t] is defined as 



(4) 



«A)) PM] = / Di()Dif)*{i(j\A\if))P[t(), t\. 



For the composition of two statistically independent sub- 
systems, which are described by two probability distri- 
butions Pi and P2 on the Hilbert spaces TLi and TL2, 
respectively, we further define the tensor product proba- 
bility distribution t] on TL = TLi <E> TL 2 as |[|] 

t] = J DipDi/j* J DtpDip* 

x V>®<p]Pi[*MP 2 [¥>,t], (5) 

where 8[-] is the functional delta distribution, and write 
as a short hand notation P[-,t] — Pi[-,t] <g) i-^hi]- The 
connection to the standard density matrix description of 
statistical ensembles is made through the identity 



p t (x,x')= / I*l>D1)*(x\il>)(il>\xr)P[il>,t]. 



(0) 



Hence, the density matrix is the covariance matrix of the 
stochastic process. 
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III. DERIVATION OF THE EQUATION OF 
MOTION FOR THE REDUCED SYSTEM'S 
PROBABILITY DISTRIBUTION 

In this section we derive the equation of motion for the 
probability distribution of the states of an open system. 
To this end, we first desc ribe b riefly the class of mod- 
els we wa nt to treat (Sec. [II A). We then introduce in 
Sec. Ill B the weak coupling assumption which enables 
us to derive an expression for the conditional transition 
probability T[ip, t + t\iPq, tp ], wher e ijj, ipo are pure states 
of the open system (Sec. IIIC2). This result will be 
further simplified by making the Markov approximation 
(Sec. [II Dj ) , and finally be u sed to obtain the desired 
equation of motion (Sec. HIE). 



Ai=A, A 2 = A\ 



A. Description of the underlying model 

Consider an open quantum system in the Hilbert space 
Til with Hamiltonian Hi interacting with an environ- 
ment consisting of harmonic oscillators, e. g., the elec- 
tromagnetic field modes. The pure states of the envi- 
ronment are elements of the Hilbert space TL 2 and the 
Hamiltonian H 2 of the environment is given by 



H 2 = },u k bj,b k , 



(7) 



where we set h = 1 for simplicity. Throughout this ar- 
ticle, bk denotes the annihilation operator for the field 
mode k = (ojk,k,\ k ) with frequency u> k , unit wave vec- 
tor k, and polarization index Afc. 

The initial condition for the state of the environment 
at time to is given by a probability distribution 



/ . 1 



dx 
2tt 



where the p a are positive, normalized weights, i. e., 
Pa > 0, = 1, 



(8) 



(9) 



and the ip a are arbitrary normalized states in TL 2 . It is 
important to note that we do not suppose that the tp a 
are orthogonal, i. e., the initial condition Eq. (J|) also in- 
cludes arbitrary mixtures of coherent states. We further 
assume that the system under consideration and the en- 
vironment are statistically independent at time to, i.e., 
that the probability distribution of the total system fac- 
torizes at time to (cf. Eq. (||)). 

The interaction of the system with its environment is 
modeled via the interaction Hamiltonian 



Hi 



5> 



Bi 



(10) 



i=l,2 



^i=-*I>fc4> B 2 = B\. 

k 



(11) 



The §k are the (real) coupling constants. For optical sys- 
tems A and A' could be, for example, the positive and 
negative frequency part of the dipole operator. It is use- 
ful to assume that the Ai are eigenoperators of Hi, i. e., 
[Hi, Ai] = LOiAi with u>i = —uj Si and oj 2 — u) s , where u> g is 
the system frequency. This is no restriction since this can 
always be achieved by choosing appropriate linear combi- 
nations of A and A* . This type of coupling describes for 
example the interaction of a two level system or a har- 
monic oscillator with a quantized electromagnetic field in 
the rotating wave approximation, which is of particular 
interest in quantum optics. Note that the restriction to 
this kind of coupling has only technical reasons. The gen- 
eralization of the presented theory to a coupling involving 
more than one system operator A, which was considered 
for example in Ref. ||, is straightforward. The Hamilto- 
nian of the total system is 



H = H 1 + H 2 + Hj. 



(12) 



In order to simplify the discussion in Sec. HID, where 



we state the conditions necessary to perform the Markov 
approximation, we will introduce at this point a decom- 
position of the Hilbert space of the environment into the 
Hilbert space 7Yb which is the state space of the bath 
and the Hilbert space TLb-? which describes the driving 
field. This decomposition reflects the fact that the elec- 
tromagnetic field which interacts with the system is in 
general produced by different sources which can have to- 
tally different effects on the system. This decomposition 
is defined as follows: 

Let /Cb be the set of modes k for which ((bk)} p 2 [- ,t ] = 0, 
i.e., the modes which do not contribute to the average 
electromagnetic field and /Cor be the set of modes k for 
which ((b k )) p 2 [- ,t ] 0. Then we define the Hilbert space 
of the bath Hb as 



H B = (g) H k , 



(13) 



where the operators Ai and Bi are defined as 



where TLk is the Hilbert space of the field mode k, and 
similarly we define the Hilbert space TLbi of the driving 
field. Obviously we have Tt 2 = TLb <£> T^Dr- We further 
assume that the states belonging to TIb and H.D r are sta- 
tistically independent, i.e., the probability distribution 
of the environment can be written as P 2 = Pb ® Pb>t, 
where Pb and Pb> t are the probability distributions on 
the bath space and on the space of the driving field re- 
spectively. This assumption is reasonable since the two 
fields are usually generated by different and independent 
sources. 

With the above definitions in mind we can decompose 
accordingly the operators Bi as follows: 
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B, = Bf 



I + I®B. 



Dr 



where 



and 



Bf = <?4> 



R B _ R Bt 



(14) 
(15) 



Dr 



B. 



Dr 



si: 



DrT 



(16) 



For the expectations of the operators Bi we thus obtain 

«%I,i] = «Oft,[.,fl 1 (17) 
i.e., the expectation of the electromagnetic field only 
depends on the state of the driving field and is inde- 
pendent of the state of the bath. For the two time 
correlation function ((B]'(s)Bj))p 2 [. i t], where Bi(s) = 
exp(i_ff 2 s)i?i exp(— iH 2 s), we find 

((Bj(s)B J )) P2[ ., ] = ((Bf \s)Bf)) PB[ ., ] 



(Sf rt ( S )i?f))p Dr[ , t] , 



(18) 

which means that this correlation functions can be split 
into the sum of two terms, where the first term only de- 
pends on the state of the bath and the second term is 
due to the driving field. 



B. Dynamics of the environment 
assumption 



weak coupling 



The basis of the weak coupling assumption is the fact 
that the environment is large compared to the system. 
Hence it is assumed that the dynamics of the environ- 
ment is given by its free time evolution; the perturbations 
due to the interaction with the system under considera- 
tion can be neglected. For the probability distribution of 
the environment this means 



p 2 [<pa =y^pa 



if ■ 



-iH 2 (t-t Q ) 



(19) 



for all t, and we have chosen the initial condition given 
in Eq. (§. 

Using the weak coupling assumption, we can further 
evaluate the expressions for the bath correlation func- 
tions ((Bf\s)Bf)) p^t] defined in Sec. |IIIA| , if the states 



of the individual bath modes are statistically indepen- 
dent, which means that the probability distribution of 
the bath is the product of probability distributions Pk 
on the Hilbert spaces Hk of the individual bath modes k. 
This will be the case for the examples we want to treat 
in this article. We obtain for example 



((Bf\s)Bf)) 



ft M] 



E 



S^ s ((&i& fc »p fc[ , to] . (20) 



It is important to note, that this correlation function only 
depends on the time argument s and the state of the bath 
at time to; it is independent of the time t. 



C. Dynamics of the reduced system 



In this section we focus on the description of the dy- 
namics of the reduced system evolving from the initial 
distribution of the total system 



P[;t Q )= P 1 [;t )® Pi\-M 



(21) 



(cf. Eq. (g) for the definition of the tensor product prob- 
ability distribution), where 



f 27r dy 

W,*o] = J ^[V'-e^o], (22) 



and P2 [ip, to] is given by Eq. (||) . Thus at time t the sys- 
tem is in the pure state ipo and the initial distribution of 
the environment is given by the probability distribution 
P%. The dynamics of the reduced system is completely 
described through the conditional transition probability 
T[ip,to + T\ipo,to], which is the probability that the sys- 
tem is in the state ip at time to + t under the condition 
that the system is in the state ipo at time to- In order 
to do this, we will start in Sec. Ill CI with the calcu- 
lation of the time evolution of the pure product state 
^a(io) = ^0 <8> <fa hi second order perturbation theory. 
Using this result we can calcul ate P[^ , t +t] and apply- 
ing a reduction formula (Sec. [II C 2] ), we finally obtain 
an exact expression for the conditional transition proba- 
bility. 



It is important to note that the conditional transition 
probability T[ip, to+r\ipo, to] also defines the dynamics of 
the reduced system, if its initial state is given by an ar- 
bitrary probability distribution Pi[if),to\. According to 
the general laws of probability theory the propagated 
reduced probability distribution Pi[ip,to + t] given the 
conditional transition probability T is 



PihMo+r] = J D^ D^T[t/j,t +T\^ ,to]Pi{^o,to]- 

(23) 



Thus it is sufficient to consider the initial conditions 
Eq. (H). 
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1. Time evolution of the total system in second order 
perturbation theory 

In this section we will study the time evolution of the 
initial distribution Eq. (|2l|) in second order perturbation 
theory. To this end, we have to calculate the time evolu- 
tion of the pure product state x f , Q (to) = ipa ® ip a - 111 the 
interaction picture which coincides with the Schrodinger 
picture at time to, the operators Ai(r) are given by 



A i (r)=A i e^\ i = l,2, 



(24) 



since A4 is an eigenoperator of H\, and the operators 
-Bj(t) are given by 



(25) 



In second order perturbation theory we obtain for the 
propagated state 

|*a(*0 + T)) = \i) ) <g> \tf a ) + ^ M^O) ® fi\Va) 

i 

+ J2 A i A ^®Sij\<Pa), (26) 



where the first and second order propagation operators 
acting on the states of the environment are 



ft = -i f ds e^'Bi{8), 
Jo 



(27) 



T pT — S 



ds / ds'e- iu ^ s+s ^e iu) ^ 
Jo 



xBjis + s'^is). (28) 

If we define the averages of fi and g^ with respect to the 
probability distribution of the environment as 

Fi = ((/i))p 2 [.,to]: Gij = {(9ij))p 2 [;t ]' ( 29 ) 
and the shifted operators fi and as 

fi = fi — Fi, 9ij = gij ~ Gij, (30) 

we can rewrite Eq. (^6|) as 

|*a(*o + r)) =L\i> ) <g> \<p a ) +^2Ai\ip Q ) <g )fi\ip a ) 

i 

+ J^A\A j \<ip )®g ij \y ot ) (31) 

i,j 

where we introduced the time evolution operator L which 
is defined as 



L = I + J2 p i A i + I] GijAlAj 



(32) 



As we will show in Sec. Ill C 2, Eq. ( |3l| ) is the basis for 
the construction of a specific reduction formula. 

We can also use Eq. ( |3l| ) to determine the time evolu- 
tion of the probability distribution of the total system if 
the initial probability distribution is given by Eq. (|2l]). 
For the propagated distribution P[^>,to + r] we simply 
obtain 



P[%to + T 



= X>«» / 
Jo 



2tt 



S^-e^Jto + r)}, (33) 



where the states ^ a (to + T ) are the propagated states 
defined in Eq. @. 



2. Reduction formula 

The key to the stochastic description of open quantum 
systems is the reduction formula ffljM. The reduction 
formula states how the conditional transition probability 
T[ip, to + t\4>o, to] is defined given the probability distri- 
bution P^i] (cf. Eq. (j33|)) on the state space of the 
total system Hi <g) H.2- In order to derive an appropri- 
ate expression for T we require the following necessary 
condition to hold: 



D^D^*(y\C <8> J|*)P[*, t + t] 
= I Di(>Di,*(il>\C\mi>,t + T\il>oM, 



(34) 

or in the more compact notation introduced in Sec. |n| 

((C ® I))p[.,t +T] = ((C))T[;t +T\A,M ■ ( 35 ) 

Note that the integral on the left-hand side of Eq. ( |34| ) ex- 
tends over the Hilbert space Tti(g)H.2, whereas the integral 
on the right-hand side only extends over the reduced state 
space Tli . The above condition ensures that the time evo- 
lution of the expectation value of an arbitrary operator 
C ®I acting on the total system's state space calculated 
using the distribution P[^f, t + r] is the same as the ex- 
pectation value of C with respect to T[ip,to + T\ipo,to\- 
However, there are infinitely many ways to construct a 
reduced probability distribution which fulfills the above 
condition, just as there are infinitely many ways to ex- 
press an arbitrary density matrix as an ensemble of pure 
states ]23| ] , and each reduced probability distribution will 
define a different stochastic process. In Ref. || it was 
shown that different processes can be formulated us- 
ing different "reduction bases" (which are bases of the 
Hilbert space of the environment). These processes were 
interpreted as the time evolution of a continuously mon- 
itored individual quantum system and each process was 
related to a specific measurement scheme p0[ . 

Here, we will follow a different approach, which has the 
advantage that it results in a particular simple, basis-free 
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reduction formula, but which is not based on a specific 
measurement scheme. If we insert the propagated prob- 
ability distribution of the total system Eq. (|33T) i n the 
left-hand side of the consistency condition Eq. (|35[) and 
neglect terms of the order g\ we obtain 



P[;to+r] 



+ Y / ^ij(M4CA j \^ ), (36) 



where we introduced the (2 x 2)-correlation matrix of the 
environment 



r 



y - ((fi fj))p 2 [-M 



= -{G ij + G* i + FZF j ), (37) 



and Fi and Gij are defined in Eq. (p9|). Note that T is 
Hermitian and non- negative. If we denote the normalized 
eigenvectors of r by in = (/in, ^ii) T , and its eigenvalues 
by Aj, and define the system (jump-) operators as 

./i = 5>w^*> (38) 

k 

then Eq. ( |36"| ) can be written as 

{(C®I)) PhtQ+T] = {M^CL\^) 

i 

= (( C ))Tl;t +r\ipoM- ( 39 ) 

Consider the following conditional transition probability: 

T[tp, t + r\ip , h] = wS[ip - w~ 1/2 Lip ] 

+ *tWiS[il> ~ w7 1/2 Jii> ], (40) 



where w — ||Z/0o|| 2 , and Wi — ||Ji-0o|| 2 - Obviously, 
the conditional transition probability defined by Eq. ( J40| ) 
satisfies the necessary condition Eq. (|35|). Since in the 
generic case the matrix T is nondegenerated we find three 
contributions to T: the system can evolve into the state 
it) -1 / 2 LiJjq, or to one of the two states w^ 1 ^ 2 JiV'o, an d 



-1/2 



J2ipo- Note, that the expression for the conditional 



transition probability Eq. ( pOj ) is exact within second or- 
der perturbation theory; however, its derivation is based 
on the special choice of the initial probability distribution 
Eq. ( pl| ) and the complete dynamics of the environment 
is contained in the functions w, Wi, Xi, and the operators 
L, and J;. 



has a large number of degrees of freedom, this formula 
is of little practical use for computations without further 
approximations. The simplest approximation scheme re- 
sulting in a conditional transition probability which is 
independent of the actual state of the environment, is 
the Markov approximation. We will use this approxima- 
tion to calculate the quantities involving the state of the 
environment, namely the ensemble average of the first 
order propagation operator Fi (cf. Eqs. ( p7| ) and (pgj)) 
and the ensemble average of the second order propaga- 
tion operato r Gg (c f. Eqs. ([2§|), and (|29|)). This will be 



done in Sec. 1IID1 and Sec. HID 2, respectively These 
results will then be combined in Sec. [II D 3 in order to 



establish an approximate expression for the conditional 
transition probability. 



1. Calculation of the ensemble average of the first order 
propagation operator 

We can simply calculate the ensemble average F{ of 
the first order propagation operator (cf. Eq. p7|)) by 
Taylor expansion. To first order in r we obtain: 



-r£t , F 2 =T£ t , 



(41) 



where £t is the mean electromagnetic field at time to, 
weighted with respect to the coupling constants gk, and 
the probability distribution P 2 [• , to] . Using Eq. ( p"7j ) we 
find 



£ ta = -i((B 



2 P 2 ];t 



= 9k((bk))p 0l [- 



.to]- 



(42) 



fcG/Cr 



Obviously, the quantities Fi only depend on the driving 
field and are independent of the state of the bath. How- 
ever, the validity of this approximation is limited to the 
validity of the Taylor expansion to first order in r, and 
this is guaranteed by the condition 



(43) 



where T£ is the correlation time of the mean electromag- 
netic field; this correlation time is of the order of mag- 
nitude of the inverse bandwidth of the spectrum of the 
driving field. 



2. Calculation of the ensemble average of the second order 
propagation operator 



D. The Markov approximation 



In Sec. pIC2 we derived an expression for the con- 
ditional transition probability T[ip,to + T\ipo,to] (cf. 
Eq. (f4(i|)) which is exact within second order perturbation 
theory. However, in situations where the environment 



Performing a Taylor expansion of Eq. (|2§|) with re- 
spect to r in order to calculate the ensemble average of 
the second order propagation operator Gij we are lead to 
a vanishing first order term, i. e., Gij is at least quadratic 
in r. However, this Taylor expansion is only valid for a 
propagation time r which is small compared to the de- 
cay time of the correlation function of the environment 



G 



((Bl(s)Bj))p 2 [.j ]. Since this correlation function is the 
sum of the correlation function of the bath and the driv- 
ing field (cf. Eq. (|l8|)), we find that Gij is the sum of 
two terms, where the first one depends only on the state 
of the bath, and the second depends only on the state of 
the driving field, i.e., 



c. — r ,a 

l J — ij 



G 



(44) 



Since we focus our interest on situations, where the decay 
time tb of the bath correlation is very small compared to 
the decay time t rj r of the driving field correlation (which 
is equal to T£ for a coherent driving field) and to the time 
scale t s of the systematic dynamics, we have to use an 
appropriate approximation of Gij for "medium" values 
of r, i. e., for values of r for which 



T B < T < TDt,T £ ,T s 



(45) 



In this case, the bath correlation function leads to a con- 
tribution Gfj linear in r whereas the driving field corre- 
lation function leads to a term Gff which is quadratic 
in t; thus we will neglect the latter. By describing the 
dynamics of the open system on a coarse grained time 
axis, i. e., by neglecting time variations which are smaller 
than tb we can perform the Markov approximation and 
extend the range of the integration over s' in Eq. (|2^) to 
infinity. The result to first order in r is 



We shall begin with the evaluation of the correlation 
matrix of the environment T: inserting Eqs. ( |4l|) and (|47] ) 
into Eq. (BT]), we find to first order in r 



r = 7T 



N-\ 
-M*e~ 



■ 1 

2iuj s to 



—Me 2iaJst ° 
N 



(48) 



By computing the eigenvalues \ and the correspon ding 
normalized eigenvectors fii of T, we obtain using Eq. (|38| ) 
the jump operators Ji. (In order to simplify the notation, 
we do not write the to-dependency of Xi and \ii explic- 
itly.) 

Furthermore, we can calculate the operator L, which 
describes the smooth time evolution of the reduced sys- 
tem, by inserting Eqs. (|l]) and @ into Eq. (||). The 
result to first order in r is 



L = I — iT^H Bl (t Q ) + H hS (t ) + H u (t ) 



(49) 



where the three operators Hm(to), Hls(^o)j an d Hu{to), 
which account for the influence of the environment on 
the open system are defined below. We can distinguish 
the effects of the two sources of the electromagnetic field: 
1. The driving field: Using Eq. (Ell) we find 



(50) 



ds 



Bf(s)B?)) PB[ . M 



(46) 



We can simplify this further bycvaluating the above in- 
tegral, as shown in Appendix |Aj. The result is: 

Gn =-r(l(N+l)-i(S + S 1 ) 

G 12 = T |Me 2l ^ to 
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G 2 i = T^M*e- 2luJsto 



Go 



(47) 



If the individual bath modes are statistically independent 
then the real constants N and Si , and the complex con- 
stant M do not depend on the time to. This would not 
be the case if some bath modes with different frequencies 
were not statistically independent. 



3. Conditional transition probability in Markov 
approximation 

In this section, we will combine the results of Sec. 



Ill C 2 , where we derived an exact expression for the time 
evolution of the conditional transitio n proba bilit y, with 
the approximations developed in Sees. Ill D 1 and III D 1 . 
to obtain the short time dynamics of the conditional tran- 
sition probability in Markov approximation. 



This is precisely the operator, which has to be added to 
the Hamiltonian of a quantum system in the presence of a 
classical electromagnetic field (according to the principle 
of minimal substitution) or a coherent driving field p4| . 
However, this result is also true for more general driv- 
ing fields, since the only assumptions we have to make 
about the nature of the driving field co ncer n its co rrela- 
tion times T£ and Trj r (cf. Sees. IIID1 and [II D 2). 
2. The bath: Using Eq. @ we find 



£ G l3 A\A 3 =: H LS (t ) + H D (t ), (51) 



where 



#ls = -(So + SJA^A + S\AJJ (52) 



is a Hermitian operator which leads to small energy shifts 
- the Lamb shift 5*0 and the Stark shift S\. The non- 
Hcrmitian operator 

Hn(t ) = -il((N + l)A*A + NAAi 

V ! \Ae 2luJ ° to - Me- 2WstQ A^A^ ) (53) 



describes the damping due to the coupling to the bath. 
If we define = Xi/jr, then the damping operator can 
also be written as 
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H o (t ) 



(54) 



Using Eq. ( p4[ ) and the fact that H^(to) is the only non- 
Hermitian operator which contributes to the operator L, 
it is clear that the weight w = \\Lipo\\ 2 is given by 



l-7r^A i <^ol4 J il^o) = l-7 



XjWj. (55) 



The latter equality ensures that to + 77- 52 • AjiUj = 1; 
which is equivalent to the normalization of the condi- 
tional transition probability. 

Summing up the above results and transforming back 
to the Schrodinger picture, we finally obtain the condi- 
tional transition probability in the Markov approxima- 
tion 



x S[ip — (I — irG (t )) M 



+ 7^^ X l \\J i j>o\\ 2 S 



4> 



(56) 



Here^we introduced the nonlinear time evolution opera- 
tor 



G(t )^ = ^eff(to)+^E A ' ; ll J ^°l| 2 ) 



(57) 



and the effective system Hamiltonian 

ffeff(to) = Hi(to) + H Dl (t ) + H LS + H n (t Q ). (58) 

Note, that the deterministic propagation with the nonlin- 
ear operator G(to) conserves the norm of the propagated 
state and thus confines the dynamics to the unit sphere 
of the Hilbert space of the system. 

The following comment may be helpful. If we combine 
the two approximations we made (second order pertur- 
bation theory and the Markov approximation) with the 
weak coupling assumption and assume that the probabil- 
ity distribution P[<3>,io] of the total system factorizes at 
time to, i- e., P[-,to] = Pi[ - >*o] ® Ah to], then the proba- 
bility distribution P[^f, t +t] of the total system at time 
to + t can be approximated by 



P[; t +T]= Pi [; t +T](g> P 2 [; t Q + T ], 



(59) 



where P\[-,to + r] is given by Eqs. ( |23| ) and (56) and 
P2[-,t + t] by Eq. (|l~9|), i.e., P[-,t + r] approximately 
factorizes. Thus we can use the same scheme for the cal- 
culation of the conditional transition probability T[ip, to+ 
t + t'\iJj, to + t] and by repeating this argument, we find 
that P[-,t] factorizes for all t. Hence, the dynamics of 
the reduced system is Markovian. 



E. Equation of motion for the reduced probability 
distribution 

One possible starting point for the derivation of the dif- 
ferential equation of motion for the reduced probability 
distribution is the infinitesimal generator Q of the pro- 
cess |25| , which is defined by its action on an arbitrary 
functional R, which may depend explicitly on time: 

g R[ ^ to] = lim TrR[M-R[M m 
t— >o r 

where the propagator V T is defined by 

V T R[i>o,t ] = J Di>Dij*R[^,to + T]T[^,to + T\4,o,to], 

(61) 

i. e., VtR^o, to] is the expectation value of the functional 
R at time to + t under the condition that the process 
starts in the state ipa at time to- Inserting Eqs. ( |5^ ) 
and (61) in Eq. ( |60|) and performing the limit r — > on 
the coarse grained time axis (i. e.,r» tb) we obtain 

J V 0% (x) 



8R[ik,ta] 



(G(t )iPo)(x) 



Ho(x) 

- [ Di;*D^(R[iP,to]-R&oM)WMi> ,t ], (62) 



where we defined the transition functional 



1 1 -MM 



(63) 



Note, that the ^"dependency of Vt / [i/ , |V , o J M is given 
through Xi and J;. Eq. ( p2| ) represents the genera- 
tor of a piecewise deterministic Markov process [ p6| , 
whose sample paths consist of deterministically prop- 
agated pieces interrupted by stochastic jumps. The 
deterministic pieces are the solution of the nonlinear 
Schrodinger equation Jt],|| 



(64) 



and the waiting time distribution function for the 
stochastic jumps is given by 



F[ipo,t ,t] = 1 - exp 



ds\\JMs)\\ 2 



(65) 



This is the probability that a jump occurs in the time 
interval [to , t) under the condition that the system is in 
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the state ipo at time to- If a jump occurs at time T then 
the state of the system at time T is 

WT)) = \\JMT-))\\ (66) 

with probability wi/ (w\-\-W2), where T_ is \m\ t ^{T — e) , 
with e > 0. 

In order to obtain the differential equation of motion 
for the reduced probability distribution, we set R[ip', t] = 
5[ip — iji'] and obtain 

V T R[iP',t]=T{<if>,t + T\i>',t], (67) 



and thus 

|- Pi ty, t] = lim - (P 1 t + r]-P 1 [if,, t] 

= hm f D ^ D r mt+T ^ t] -^-^ ] m',t] 

T-»0./ T 



(68) 



If we insert Eq. ( |62| ) in Eq. (|6q ) and evaluate the inte- 
gral we obtain the differential equation of motion of the 
reduced probability distribution 



(G(t)V>)(a 



+ I Di//Di//*{w[il>\il>', t]Pi W, t] - WWW, t\Pi[^ t] } 



(69) 

which is a Liouville-Master Equation . Eq. (^9|) implies that for the conditional transition probability the integro- 
differcntial equation 

& f f S 



(G(t)i>)(x)- $^(G(W) {x)\T[i>,t\^ M 



(70) 



holds for any t > to- This equation can be solved either 
analytically for some simple systems using the Hilbert- 
space path integral representation of the process [|l2|j27| 
or numerically by either simulating the_process [|7l|8]l or 
using a recursive computation scheme 



IV. REDUCED HEISENBERG PICTURE 
OPERATORS 

Studying the dynamics of open systems, one is not 
only interested in the time evolution of expectation val- 
ues of some system operator X but, moreover, one wants 
to investigate the time evolution of an arbitrary matrix 
element of this operator and obtain quantities such as 
((f>o,tQ\X(T)\ipo,t ). Note that for an open system this 
is only a short hand notation. What we really mean by 
this is 



b ,to\X(T)\i> ,t ) 

^5> Q ($ Q ,to|e 1 ^ (X 



7)e 



l*a,*0>, (71) 



where H is the Hamiltonian of the total system (see 
Eq. (12)) and the states $ Q and $„ are defined as 
&a — 4>o® l Pa and \& Q — ^o®Lp a i i. e., we assume, that the 
state of the environment is given by the probability dis- 
tribution Eq. (||). We will refer to this quantity (Eq. (fn])) 
as the matrix element of the reduced Heisenberg picture 



operator X(t) and use it in Sec. V B 2 for the calculation 
of time-ordered multitime correlation functions. 

In the case of a closed system, (<fr , t \X(T)\tp , to) 
is easily calculated by solving the corresponding 
Schrodinger equation with the two initial conditions 
(j>(to) — (f>o and i/)(to) = V'o separately to obtain <fr(to + t) 
and ift(to + t) and then evaluating the scalar product 
(<f>(to + r)\X\ip(to + r)). A naive generalization of this 
scheme to an open system - which will lead to wrong 
results - would be the following: instead of solving the 
Schrodinger equat ion, simulate the stochastic process de- 
fined in Sec. 1IIE| with the two initial conditions <j)o and 
ipo to obtain <j>(to + t) and i/>(to + t). Then evaluate 
the scalar product (<^>(to + r)\X\ip(to + r)) and average 
over a sufficiently large ensemble of realizations. How- 
ever, what we really obtain using this procedure is the 
quantity 



x T[&to + T|<fo,to]T[^to + T|^oM (72) 

which is, in general, not equal to (</> , to\X(T)\ipo, to)- 
This can be most easily seen by considering the spe- 
cial case 4>o — tpo'- According to the defin ition o f 
the conditional transition probability (cf. Sec. [II C 2] ), 
T[ip,to +T\ipo,to] has to fulfill the necessary condition 



(V>o,to|*(r)|Vo,io) 



Dif>Dip*(<ip\X\ 







x T[il),t + T\ifa,to]. 



Obviously, the quantity f(tpo,ipo,to, T ) defined in 
Eq. ((72)) is, in general, different from the right-hand side 
of Eq. (|73|), and thus this approach will lead to wrong 
results. 

However, as we will show below, it is possible to 
construct a stochastic process which describes the time 
evolution of 4>o and ipo simultaneously and whose sam- 
ple paths can be used to estimate the desired quantity 
(0o, to\X(T)\ip , t ). In order to do so, we describe the 
dynamics of the open system in the doubled Hilbcrt 
space Hi = Hi © Hi, such that the "state" of the open 
system is given by a normalized pair of state vectors 
8 = c- 1 / 2 (0,V) T , where c = \\<j>\\ 2 + \\tp\\ 2 and T de- 
notes the transposed vector. Accordingly, we define the 
extension of a system operator X to the doubled Hilbert 
space Hi as 



X 



I 
X 



Using the raising operator 

S+ = 



/ 




(74) 



(75) 



whose expectation value c(8\S + \8) in the state 8 is the 
scalar product (<f>\ip), we can write the matrix elements 
of X as 



(<p\X\i/>) = c(0\X*S + X\0). 



(76) 



This equation is the crucial point of the following con- 
struction: it allows to write matrix elements of re- 
duced system operators as expectation values in the dou- 
bled Hilbert space Hi. The same idea will be used in 
Sec. V B 2 for the development of an algorithm that 



enables the determination of arbitrary correlation func- 
tions. Note, that the scalar products on the left-hand side 
and on the right-hand side of the above equations are de- 
fined in different Hilbert spaces, namely in Hi and Hi, 
respectively. (The on the right-hand side of Eq. (frfj) 
is introduced to hi ghli ght the connection of the above 
equation with Eq. (105) (see below)). We also introduce 
a reduced probability distribution Pi [8, t] on the set of 
normalized states of Hi , which is the probability density 
of finding the system in the state 8, and a conditional 
transition probability T[8, to + t\8q, to] which is the prob- 
ability density that the system is in the state 8 at time 
to + t under the condition that the system is in the state 
8q at time to. 

The goal of this section is to derive in analogy to the 
scheme we presented in Sec. Ill the equation of motion 



for the conditional transition probability T, which has to 
fulfill the necessary condition 

(0o, *o|*(r)HMo) = c o((^S+X))~ lto+Tleoto] , (77) 

where c = ||0o|| 2 + \\ipo\\ 2 and 8 = c o ~ 1/2 (0 o , ipo) T ■ 
The above condition states that the matrix element 



(73) (0o, to|-X"(r)|^o, *o) 01 a reduced Heisenberg picture oper- 
ator is proportional to the expectation value of the oper- 
ator X'S + X with respect to the probability distribution 
T . In order to achieve this we have to make the follo wing 
assumptions, which are similar to those made in Sec. Ill: 
i) According to the weak coupling assumption, we will as- 
sume that the time evolution of the environment is given 
by its free evolution (cf. Eq. Jl9|)). ii) The probabil- 
ity distribution of the total system at time to, P[0,to] 
(which is defined on the Hilbert space H = Hi ® H%), 
factorizes and is given by 



p[e,t ] 



Pa 



dx 
2n 



6[G 



% ® fa 



(78) 



iii) Second order perturbation theory and finally iv) the 
Markov approximation. If we combine these approxima- 
tions we find after some calculations which arc similar to 



those presented in Sec. Ill, that the necessary condition 
Eq. ( [77[ ) to first order in r (and to second order in the 
coupling gk) reads: 

(0o,io|*(r)|Vo,io) = (M^XL^o) 

i 

±Co%&S+x))f y (79) 



where L, and A;, and J.; are defined as in Sec. HID 2 



Note that this equation is completely analogous to 
Eq. (p9|). This leads us directly to a conditional tran- 
sition probability to first order in r 



T[8, t + t\9 , t ] = w5[8 - w^^Ldn] 



IV, 



Wi 1/2 Ji8 ], (80) 



where we have defined 







7r f + \ _ ( H e g(t) 
Hes{t) ~ { H cff (t) 



J t 
J; 



(81) 



The effective Hamiltonian H c g(t) is defined in Eq. (pq), 
and L = I — irH c g. The appropriate normalization con- 



stants are w — ||L#ol| 2 an( i Wi = II J, 



Eq. 



we find 



In analogy to 



w 



7T 



(82) 



which ensures the normalization of the conditional tran- 
sition probability T. 

The equation of motion for the conditional transition 
probability T[9, t\6o, to], where t > to, ca n be found us- 
ing the same scheme as described in Sec. HIE: Calcula- 
tion of the generator Q and its action on the functional 
R[8',t] = 5{6 — 8'] leads to the Liouville-Master equation 
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(G(t)e)(a 



69* (x) 



G(t)6) (x)\T[0,t\0 o ,t o 



DQ 



'D9'*{w[e\9', t]f[6', t\9 , t ] - W[9'\9, t]T[9, t\6 , to]}, 



(83) 



where we defined the functional Wirtinger derivative with 
respect to states belonging to Tti as 



/ 



89(x) 



8<f>(x) 
S 

Km*) J 



(84) 



The generator for the continuous time evolution is de- 
fined in analogy to Eq. (pl\) as 



G(t)9= (tfeffW+z^A^lJ^H 2 ^ 9, 



(85) 



and the transition functional is given by (cf. Eq. (p3)) 



w[e\e ,t] =iJ2 x > 



-1/2 



Ji@0 



(86) 



Again, the stochastic process which is defined through 
Eq. ( |83| ) is a piecewise deterministic Markov process, 
where the deterministic pieces are the solution of the non- 
linear Schrodinger equation 



ij t 9 = G(t)9 



(87) 



and the waiting time distribution for the stochastic jumps 
under the condition that the system is in the state 9q at 
time to is given by 

F[9 Q ,t Q ,i\ = l-exp(- 7 Y^\ i £ dsWJ^W^j . (88) 



The Liouville-Master equation fl83| ) also determines the 
equation of motion for the reduced Heisenberg picture 
operator X(t): Expanding X(t) in terms of its matrix el- 
ements and inserting the Liouville-Master equation (^) 
we obtain after some calculations 



j f X(t) = i[Ht + H LS + Hur, X](t) + ~ 

i 

x {2 (j}XJ>) (t) - (j}JiX) (t) - (XJ}J % ) {£)}. (89) 

Thus the equation of motion of the reduced Heisenberg 
picture operator is the adjoint equation of the quantum 
master equation pO|,pl| and we find 



^XVit, t ){ Po }} = Tr sys {x(t){p }} (90) 



for any initial density matrix po, where V(t, to) is the time 
evolution super operator of the quantum master equa- 
tion. This relation is the basis of the quantum regression 
theorem |J|Ji]J. 

We close this section with a few remarks: For the ini- 
tial condition #0 = (V'O) ipo) T /V^ the stochastic process 
defined in t his wa y reduces to the stochastic process de- 
fined in Sec. Ill E ; formally this can be expressed through 
the identity 



T[9,t\9o,t }=S 



where T\4>, tUpp, to] is the solution of the Liouville-Master 
equation (|69|). Thus, the dynamics of the system de- 
scribed in the doubled Hilbert space Hi reduces to 
our earlier description of the dynamics in the ordinary 
Hilbert space "H\. This is an important property of the 
conditional transition probability T which establishes the 
link between the two approaches. 

For the numerical simulation of the stochastic process, 
it is useful to omit the normalization requirement and 
use the linear operator H c r instead of G as the genera- 
tor for the deterministic motion. Thus, for the numer- 
ical simulation we use an unnormalized state vector 8, 
and the deterministic pieces are the solution of the linear 
Schrodinger-type equation 



dt 



— Hcfft 



(92) 



and the waiting time distribution for the stochastic jumps 
for unnormalized states is given by 



F[8o,to,t] = l-\\§(t)\\ 2 , 



(93) 



under the condition that the state 60 is normalized and 
6(t) is the solution of the above Schrodinger-type equa- 
tion with the initial condition 9 (to) = 9q. This equation 
is easily proven by differentiating Eqs. ( |88| ) and (|9^) with 
respect to t and comparing the results. The above pro- 
cedure has the major advantage, that the Schrodinger 
equation (|9^) for 9 = (0, -0) is linear and that the two 
components are decoupled. This leads to a large reduc- 
tion of the time necessary for the calculation of the de- 
terministic time evolution. In addition, we do not have 
to calculate the waiting time distribution separately. 

In Ref. (2^] Gisin formulated the Heisenberg picture 
of the quantum state diffusion model of open systems. 
In analogy to our approach, he used a coupled pair of 
stochastic differential equations for the calculation of 
Heisenberg picture operators. His approach has the ad- 
vantage that the scalar product of two state vectors <j> 
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and ip is constant in time, i. e., the matrix elements of the 
identity operator are calculated correctly in each realiza- 
tion of the stochastic process, whereas in our approach 
the identity operator is treated like any other system op- 
erator, i.e., for the calculation of the matrix elements 
one has to average over many realizations. On the other 
hand, Gisin's approach is limited to the calculation of 
matrix elements with respect to nonorthogonal state vec- 
tors, i. e., the stochastic differential equations are not de- 
fined if {<j)\i}}) — 0. Furthermore, the quasi linear stochas- 
tic equations proposed for the numerical simulation, are 
not decoupled - in contrast to our approach - which is a 
disadvantage from a numerical point of view. 



V. CALCULATION OF MULTITIME 
CORRELATION FUNCTIONS 



A. Symmetric time-ordered multitime correlation 
functions 



The theory we presented in Sec. Ill allows the cal 



culation of one time expectation values of an observable 
such as (ipo, t \X(T)\ip , t ). This theory can easily be 
extended to the calculation of symmetric, time-ordered 
multitime correlation functions such as 



expression for g(ipo,to,ti,t 2 ) in terms of probability dis- 
tributions defined on the Hilbert space of the total system 
can be written as 



g(ip ,t ,tx,h) = 

X (*i|X t X|*i)f 



DfiD^l J D^> 2 D%(^ 2 \Y\^ 2 ) 



P[*i,*i]. (96) 



The time-dependent probability distribution P and the 
conditional transition probability T for the total system 
can be calculated using the unitary time evolution, 

P[*i,ti] = J2p a S[^i - exp(-iH{t - ti))ipo <g> ip a ] 

Oi 

f , [* 2 ,t2|*i,ti]=*[*2-exp(-tir(t2-ti))*i], (97) 

where H is the Hamilton of the total system (cf. Eq (fL2|)). 
If we apply the scheme we presented in Sec. |ll (i.e., 
the weak coupling assumption, second order perturbation 
theory, the reduction formula, and the Markov approxi- 
mation) to express P and T in terms of the conditional 
transition probability T of the reduced system we find 
that g(ipo, to, t\, t 2 ) is approximately given by 



D^DiPl / D^ 2 DrMY\ip2 



g(i>o,t ,ti,t 2 ) = (VMo|* + (ti)y(i 2 )X(ii)|^,io), (94) x \\X^\\ 2 T 



Will' 



ThMilVMo], (98) 



where X and Y are arbitrary operators and to < t\ < t 2 . 
As in Sec. IV we use the short hand notation 



(Vo,to|A't(t 1 )y(t 2 )X(t 1 )^o,to) 

= Y,p^^Mx ] {ti)y{^)x{t x )\^ a ,t ), (95) 

a 

where ty a — ?Ao <8 (fi a is a pure product state of the total 
system. Thus, at time to the system is in the pure state 
"00 and the state of the environment is given through 
the probability distribution P 2 (cf. Eq. (pf)). An exact 



where the ipi are states of the system under consideration 
and the conditional transition probability T satisfies the 
Liouville-Master Equation Eq. (|7(i|). For the stochastic 
process which is used to simulate Eq. ( pl| ) this means the 
following: the process starts in the state %po at time to 
and is propagated to the state ip\ at tim e t\ us ing the 
stochastic time evolution defined in Sec. IIIC2. Then 



the process jumps to the normalized state X^i/||X-(/>i|| 
and is propagated to the state ip 2 at time t 2 . 

For the most general symmetric time-ordered correla- 
tion function we obtain in a similar way 



g(i>o,t , ...,t n ,t n+ i) = (^oMA{h)---Xl(t n )Y{t n+1 )X n {t n )---X 1 (t 1 )\^oM) 
= J ' DipiDtpl ■ ■ ■ J D^ n+ iD%j;* n+1 {'tp n+ i,t n+ i\Y\-ilj n+ i,t n+ i)wi ■ ■ ■ w n x 

x T[t[) n+ i,t n+1 \w n ~ 1/2 X n ip n ,t n ] ■ ■ ■T[^2,t2\wi~ 1/2 Xiijj 1 ,t 1 \T[il; 1 ,t 1 \ip ,t \, 



(99) 



where Xi and Y are arbitrary system operators, and 
Wi = \\Xiipi\\ 2 , and to < t\ < ■ ■ ■ < t n +\- This class 
of correlation functions contains for example the cor- 
relation (a + (t)(j + (t + t)o~ (t + t)ct~(<)) which has been 
measured in the famous Hanbury-Brown and Twiss ex- 
periment (see, e. g. [p8[) or, more recently, for a single 
ion in a trap by Diedrich and Walther 129] . 



The stochastic simulation algorithm which we use to 
compute symmetric time-ordered multitime correlation 
functions can thus be summarized as follows: 1. Start 
in the state ipo at time to and propagate up to the 
state ipi at time t\ using the stochastic time evolution 
(cf. Sec. |III Ej ). 2. Jump to the state AiVVII^i^ill- 
3. Propagate the state Ai?/;i/||Xi-0i|| at time t\ up 
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to the time t2, etc. 4. Compute the expectation value 
(<ip n+ i,t n+ i\Y\ip n+ i,t n+ i). By repeating this procedure 
sufficiently often, we can estimate the unknown correla- 
tion function by averaging over the sample paths. 

In the preceding discussion we limited ourselves to a 
pure initial state ipo of the system. However, this is not a 
real restriction: if the initial state of the open system is 
described by a probability distribution Pi[ip,to], we find 

«Xj(ti) • ■■Xl(t n )Y(t n+1 )X n (t n ) ■ ■■X 1 (t 1 ))) Pl[ . M 
= f Dil>D4*P l ty,t a ]g(!l>,t a ,...,t ni t n+1 ) (100) 

Thus, we can simulate this type of correlation functions 
by drawing a random initial state tp with probability 
PiH'yto] an d using the above simulation algorithm. 



B. General time-ordered multitime correlation 
functions 

Especially in the context of quantum optics, measure- 
ments are not restricted to obtain information about 
symmetric time-ordered correlation functions. For ex- 
ample, the fluorescence spectrum of a two level atom is 
the Fourier transform of the stationary expectation value 
(((t + (t)o-~)) s which is a special case of the more general 
time-ordered multitime correlation function 

9 (ipo, h, h, —,t n , si, s m ) 

= (xPoMX^h) ■ ■ ■ Xl(t n )Y m (s m ) ■ ■ •yi(s 1 )|Vo,*o>, 

(101) 



i 

X^MY=-Ye- 2 ^ N 

fe=i 

(x + e 27Tk ^ N Y^ M (X + e 2 * kt / N Y^ , N > 3 (102) 

(which yields the polarization identity |l9| for N = 4). 
This identity is easily proven by computing the product 
and recognizing that the prefactor of every term in the 
sum is a geometric series. For example the correlation 
function ((X^ (t + r)Y(t))) could be rewritten as 



((Xl(t + T)Y(t))) = -Y,e- 2 " kl/4 
fe=l 

(( (l + e 27rW / 4 y(t)V X\t + r) (l + e^ kl/i Y(t)) V>. 

(103) 

This is a sym metric time-ordered correlation function. 
Relation ( |l03| ) was used by Dalibard et al. jf) for the com- 
putation of the spectrum of fluorescence of a two level 
atom. However, for the calculation of a general time- 
ordered correlation function, one would need to calculate 
at least 3 n + m -' r - s - 1 symmetric time-ordered correlation 
functions, where r is the number of indices i for which 
tj is equal to some Sj and s is the number of indices i 
for which ti and Xi are equal to some Sj and Yj , respec- 
tively. Especially for higher order correlation functions, 
this procedure is very time consuming. 



where to < • • • < t n , and to < si < • • • < s m , and Xi and 
Yi are arbitrary system operators. In fact, according to 
Gardiner jl9| and Gardiner and Collett this kind of 
multitime correlation function is the only one that arises 
in the quantum theory of measurement. 

For the calculation of this time-ordered correlation 
function we will present two distinct approaches: we 
can express the general time-ordered correlation func- 
tion through a linear combination of sy mmetric time- 
ordered correlation functions (Sec. VB 1 ) or we can de- 
fine a stochastic process in the doubled Hilbert space 
Tii © Hi (cf. Sec. IV) which allow s a direct calculation 
of the sought quantity (Sec. VB 2] ). 



1. Reduction to symmetric time-ordered correlation 
functions in Tii 

In order to express the gene ral time-ordered correlation 
function defined in Eq. (101) as a linear combination of 
symmetric time-ordered correlation functions as defined 
in Eq. (M), we can use the operator identity 



2. Calculation by the stochastic process in the doubled 
Hilbert space Tii 



In Sec. V A , we demonstrated how the stochastic pro- 
cess which we use to simu late o ne-time expectations of 
system operators (cf. Sec. [II E ) could be applied to the 
simulation of symmetric time-ordered multitime correla- 
tion functions. Similarly, the stochastic process which 
we use to calculate matrix elements of a reduced Heisen- 



berg picture operator (4>o, to\X(T)\ipo, to) (cf. Sec. |IV|) 
provides a direct method for the calculation of general 
time-ordered correlation functions of the type given by 
Eq. (101). The reason for this fact is — as we will show 
below — that these correlation functions can be expressed 
by symmetric correlation functions in the doubled Hilbert 
space Hi in the following manner: Order the set of 
times {t±, ■ ■ ■ t n , Si, • • • s m } and rename them such that 
7*1 < • • • < r n+m -k where k is the number of coinciding 
time points. Then define a set of Schrodinger operators 
Fi on Hi as 
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x t 

I 



if n =ti 7^ Sj for some i and all j, 



F, = < 



I 




Y \ , if ri = sj 7^ for some j and all i, 



X 

g 1 F J , if n = ti = sj for some i and j. 



(104) 

Using this definitions and the definition of the raising 
operator S + (cf. Eq. (|75|)) we find from Eq. ( 101 ) 



g(ipo,to,ti,...,t n ,si,...,s m ) 

= c Q (e a \F}( ri ) ■ ■■Fl{r q )S + F q {r q ) ■ ■ ■ F x {r x )\0 o ), 

(105) 



where q = n + m 
C^ 1/2 ^0,^0). Eq. £oj| 

multitime correlation function 



fc, c = 2||^ || 2 , and O = 

is a symmetric time-ordered 
that can be calculated 



applying the simulation algorithm in a doubled Hilbert 
space described in Sec. VA. Note that for the special 
choice n = 0, m = 1, t\ = to, Co = \\4>\\ 2 + \ \ip\\ 2 , and 
9 a = c o " 1 ^ 2 (0, tp) the right-hand side of Eq. ( |105| ) is just 
the matrix element ((f>\Y\i(}) of some system operator Y 
(cf. Eq. ©). _ 

As an explicit example of the general formalism 
we will consider the two time correlation function 
{ipo,to\X^(ti)Y{t 2 )\^o,to), with t < h < t 2 : for this we 
obtain the following simulation algorithm: 1. Start in the 
normalized state ipo at time to and propagate to the state 
ipi at ti me t\ using the stochastic time evolution in 7ii 
(cf. Sec. Ill E| ). 2. Jump to the normalized pair of state 
vectors c 1 " 1/2 (X^i, ipi), where c\ — ||X^i|| 2 + ||^i|| 2 . 
3. Propagate to the state (4>2,4>2) at time ti using the 



'(Xipi,^i) 
Propagate to the state (02, "02) 
stochastic time evolution defined in the doubled Hilbert 
space (cf. Sec. IV). 4. Evaluate the scalar product 
(4>2\Y\i/j2) and weight it with the factor c%. Repeat this 
procedure sufficiently often and average over the different 
realizations. Using Eq. (^) it is straightforward to show 
that applying this scheme we determine the quantity 



the calculation of the fluorescence spectrum of a driven 
two level atom as an example. 



VI. EXAMPLES 

A. Coherently driven two level atom 

Consider a two level atom with the Hamiltonian Hi = 
uj s (T + a~ , where are the pseudo spin operators of the 
atom, driven by a coherent field. The state of the envi- 
ronment is given by 



|^)=n^(^e--*)|0) 



(107) 



with probability 1, where D is the unitary displacement 
operator and f3k is the amplitude of the coherent excita- 
tion of the field mode k. We will assume that the band- 
width of the coherent excitation is small, i.e., Eq. ( |4^ ) 
holds. Using Eq. (42) we find for the mean of the elec- 
tromagnetic field 



St = -%((B 



2/)p 2 



-iu k t 



(108) 



Bl- 



and the driving term in the effective Hamiltonian is thus 
given by 

#Dr= ]T g k (p k e-^ t a+ +Pia-e^" t ). (109) 

Since the bath is in the vacuum state, we find using 
Eqs. (|A|) and ( |Al2j ) from the Appendix 



n(u> s ) — m(uj s ) = 0. 



(110) 



Computing the eigenvalues of the correlation matrix of 
the environment T (cf. Eq. (|48"|)) leads to Ai = 1, and 
A2 = 0, and thus we find the jump operators J\ = o~ , 
and J2 = c + . Neglecting the Lamb shift we obtain for 
the effective Hamiltonian 



Tr 



sys 



rU(t 2 ,ii){u(ii,M{IV'o)(V'o|}X t }}, (106) H ef£ -(uj s -i^)a + a 



where V(t',t) is the propagator of the quantum mas- 
ter equation. This is in complete agreement with the 
standard definition of the two time correlation function 
(ipo,to\X^(ti)Y(t2)\ipo,to) and, hence, in complete 
agreement with the quantum regression theorem. 

The above algorithm is easily generalized to the calcu- 
lation of higher order correlation functions. Note that by 
using this approach we have to calculate only one sym- 
metric time-ordered multitime correlation function in a 
doubled Hilbert space, instead of 3 n +™- r - s - 1 correlation 
functions that have to be calculated using the method 
described in Sec. VB 1. We will compare the numerical 



+ gk(Pke-^ kt a++P* k a-e l ^ t ). (Ill) 

This is the same effective Hamiltonian, that is used in 
other stochastic wave function approaches for the coher- 
ent time evolution of a system coupled to the vacuum, 
see, e. g. Refs. |,|||. 

We can now use the simulation algorithm described 
in Sec. VA to calculate the correlation function 



performance of both methods in Sec. 



VIA where we use 



((er+ (t)a + (t+T)a~ (t + r)a~ (t))) s which is the probability 
of emitting a photon at time t and a second photon at 
time t + r in the steady state. This correlation function 
could be measured in a Hanbury-Brown and Twiss ex- 
periment. We do this by drawing a random initial state 
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from an uniform distribution on 7ii and then propagat- 
ing this state for a time t — 307 _1 in order to reach the 
steady state regime. Then we apply the above described 
simulation algorithm to calculate the sought quantity. In 
Fig. ^ we compare the numerical results with the analyt- 
ical solution obtained by the quantum regression theo- 
rem (solid line) for a resonant coherent excitation with 
St = 5je~ luJst (which yields the Rabi frequency ft = IO7) 
and 10 5 realizations. The error bars indicate the esti- 
mated error calculated at each point using the estimated 
standard deviation of the sample of realizations. Obvi- 
ously, our calculations are in good agreement with the 
analytical results. 

We can also compute the spectrum of the fluorescence 
radiation which is the Fourier transform of the correla- 
tion function (((J + (t + T)a~ (t))) u sing the simulation al- 



together with Eq. ([42]) we find for the mean of the elec- 
tromagnetic field 



gorithms presente d in S ec. VB 1 (denoted by Method I 
in Fig. I) and Sec. |VB2| (denoted by Method II in Fig. |) 
respectively. For the same parameters as above, both al- 
gorithms converge to the analytical solution, first given 
by Mollow pi]] . However, the numerical performance of 
both algorithms is quite different: Fig. || shows the nec- 
essary CPU time on an RS6000 workstation in order to 
achieve a desired numerical precision in a log-log plot. 
The solid lines represent the mean square deviation of 
the numerical solution from the exact solution and the 
dashed lines show the mean estimated standard deviation 
of the numerical solution. Obviously the latter quantity 
provides for both algorithms a very good measure of the 
accuracy of the numerical simulation. A closer analysis 
reveals, that Method II is for a given accuracy by a factor 
of 3 faster than Method I. Actually this result is better 
than we expected, since the time Method I uses for the 
computation of a single realization is only by a factor of 
1.7 greater than the time Method II needs. Thus, the ac- 
curacy of a single realization is higher for Method II than 
for Method I. We emphasize that, in order to obtain a 
realistic picture of the scaling of the CPU times, we used 
a numerical integration of the corresponding differential 
equations (Eqs. ( J64| ) and (|87|)), instead of the well-known 
analytical solution. 



B. Two level atom in a squeezed vacuum 

As a second example, we want to investigate a two 
level atom with the Hamiltonian Tii = Lu s a + a~ , driven 
by a displaced squeezed vacuum. Thus, the state of the 
environment is given by 



W) = J{D{p k e-^ t )S(T k e- ii ^)\Q) 



(112) 



with probability 1, where D is the unitary displacement 
operator and S is the unitary squeezing operator. Using 
the identities pl| 

S\re~ 2l ' t >)b k S{re- 2l ' l >) = b k coshr - 6£e -2i * sinhr 
S^re-^blSire- 2 ^) = b{ coshr - 6 fc e 2 ^ sinhr (113) 



St 



-i((B 



■ £ 

fce/c Dr 



gk(3 k e 



-iujfet 



(114) 



and thus for the driving term in the effective Hamiltonian 
H Dl = ]T g k (p k e- iu " t a + + p* k a-e iu " t ). (115) 

It is important to note that the driving term is completely 
unaffected by the squeezing and again the coherent ex- 
citation acts like a classical field. If we assume that the 
bandwidth of the coherent excitation is small and that 
the squeezing is homogeneous over a wide bandwidth 
near the resonance frequency uo s then Eq. (^) holds and 
we can calculate the parameters N and M a ppear ing i n 
the bath correlation matrix V using Eqs. (A2) and (|A12|) . 
We find 

7V = esinh 2 r s , M = ee~ 2i0s sinhr s coshr s , (116) 

where r s and <p s are the parameters of the squeezing at 
the resonance frequency uj s and the efficiency parameter 
e is defined as e = | cr| / 4vr where a is the solid angle in 
which the squeezing occurs. If the squeezing is homoge- 
neous over the complete solid angle An, then e — 1 and we 
find a perfect squeezed vacuum with \M\ 2 = N(N + 1) 
(minimum uncertainty state); if only a finite solid an- 
gle is squeezed, then we find \M\ 2 = N(N + e), which 
corresponds to an imperfect squeezing. Thus a perfect 
squeezing over a finite solid angle leads to the same dy- 
namics as an imperfect squeezing within the complete 
solid angle An. 

In order to obtain the rates jXi and the corresponding 
jump operators Ji which are necessary for a stochastic 
description of the system we have to calculate the eigen- 
values and eigenvectors of the bath correlation matrix T. 
In terms of N and e we find 



Ai, 2 = N + \ ± y/ N(N + 7) + I 

Ji = co&ee^'A- smOe-^'Ai 
J 2 = sin (9e i0s A + cos 6»e^ i0s A f , 



(117) 



where tan26> = 2y/N(N + e). 

The spectrum of fluorescence can be obtained by calcu- 
lating the correlation function ((ct + (t)(t _ )) s in the steady 
state, and then taking the Fourier transform of this quan- 
tity. In Fig. |§| we present numerically calculated spectra 
of resonance fluorescence for a coherent driving field with 
a Rabi frequency of ft = IO7 (S t = ^ ie -^^-<P^)) and 
various squeezing parameters (solid lines) and compare it 
with the vacuum spectrum obtained by Mollow pl| (thin 
line). In Ref. |52| it was shown that in the strong field 
limit the spectrum of resonance fluorescence is sensitive 
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to the relative phase of the driving field and the phase of 
the squeezing, i.e., it depends on <j> = 2{<f> s — 4>l): For a 
perfect squeezing and cf> = the linewidth of the central 
peak is reduced to a subnatural level, whereas the side- 
band peaks are broadened. This is shown in Fig. [jj (a) 
and H (c). In the limit e — > but TV = const, (i. e., very 
strong squeezing of only a few modes) the shape of the 
central peak is approximately identical with the shape 
of the central peak in the Mollow triplet; the sideband 
peaks are still broadened (cf. Fig. || (b) and || (d)). On 
the other hand, for <f> = ir the central peak and the side- 
band peaks are broadened for both e = 1 (cf. Fig. || (e) 
and | (g)) and e = (cf. Fig. | (f) and | (h)). 

For the numerical calculation of the correlation func- 
tion (((j + (t)(7~)) s we used the stochastic simulation 
algor ithm in the doubled Hilbert space described in 
Sec. VB 2| and 10 4 realizations. Note, that it is use- 
ful from a numerical point of view to subtract the 
limit lim T _ >00 ((cr + (T)(T _ )) s from the calculated correla- 
tion function before taking the Fourier transform, since 
this constant corresponds to the 6— shaped coherent part 
of the spectrum, which is neglectable in the limit of 
strong driving fields. The agreement of the numerical 
with the analytical solution given in Ref. |32| is very good 
for all the examples presented above (the relative error 
of the correlation function is of the order 10 -2 ). Note 
that the analytical solution was obtained by assuming 
that all modes of the electromagnetic field are squeezed. 
However, the above discussion makes clear that it can 
also be used for the more general case considered here by 
setting \M\ 2 — N(N + e), i.e., by assuming an imperfect 
squeezing. 



C. Thermal mixture of coherent states 

As a last example we want to investigate an open sys- 
tem coupled to a heat bath in a thermal state, which is 
represented by a mixture of coherent states. Thus the 
probability density Pb [<p } to] of the bath is given by 



PbIM = (g) ft Mo], 

keic B 



(118) 



where 



(-M 2 ( e ^ /fcT -lJjSfy- e lx \a)} , (119) 



x exp i 



and | a) is an eigenstate of the destruction operator 
with com plex e igen value a. By inserting Eq. (HE) in 
Eqs. (A2) and (|A12|), respectively, we obtain 



N = 



1 



o k /kT _ y 



M = 0, 



(120) 



as expected, and hence the jump operators Ji = Ai and 
rates 7A1 = j(N+l) and 7A2 = jN. Note that we would 



have obtained the same values for TV and M and hence 
the same equation of motion for the reduced probability 
distribution by representing the probability density Pb 
of the heat bath in terms of eigenstates of the number 
operator The reason for this invariance is that all 

parameters which describe the influence of the environ- 
ment on the open system are expectation values of linear 
environment operators and hence independent of the ac- 
tual representation of the mixture. 



VII. SUMMARY 



The primary goal of the present paper was to show 
how stochastic wave function methods have to be gener- 
alized in order to allow the systematic evaluation of the 
matrix elements of reduced Heisenberg operators and of 
arbitrary multitimc correlation functions. The starting 
point for the derivation of the equation of motion for the 
reduced system's probability distribution is the unitary 
time evolution of the total system (i.e., the system to- 
gether with its environment). From this, we infer the 
dynamics of the reduced system through the following 
necessary condition: the time evolution of the expecta- 
tion value of an arbitrary system observable C calculated 
using the reduced probability distribution has to be equal 
to the time evolution of the expectation value of C ® I 
in the total system (up to second order perturbation the- 
ory). By making the week coupling assumption and the 
Markov approximation we can further simplify the ex- 
pression for the reduced system's probability distribution 
and finally obtain the desired equation of motion. This 
equation is a Liouville-Master equation and hence the 
underlying stochastic process is a piecewise determinis- 
tic Markov process. We emphasize, that the only major 
assumptions we have to make about the state of the en- 
vironment concern the relation of the time scale of the 
reduced system's dynamics compared to the time scale 
of the dynamics of the environment. 

The basic idea underlying the approach to the evalu- 
ation of multitime correlations presented here is the in- 
troduction of a doubled Hilbert space. The latter allows 
to write matrix elements of reduced Heisenberg picture 
operators as expectation values. Employing this ansatz 
a stochastic process in the doubled Hilbert space is eas- 
ily constructed which directly leads to the determination 
of arbitrary correlation functions. The process in the 
doubled Hilbert space is again a piecewise deterministic 
process described by a Liouville Master equation for a 
probability distribution on the doubled Hilbert space. 

To complete this article we have discussed three typi- 
cal examples of quantum optics: (i) A coherently driven 
two level system coupled to the vacuum. This example 
is mainly used as a basis for a numerical comparison of 
two simulation algorithms for the calculation of multi- 
time correlation functions, namely our approach and an 
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algorithm proposed in Ref. EJ. For this example our ap- 
proach was faster by a factor of 3. (ii) A coherently driven 
two level system coupled to a squeezed vacuum. In this 
example we derive the stochastic time evolution of a sys- 
tem coupled to a squeezed vacuum which extends over a 
finite solid angle (i.e., not over the complete solid angle 
47r) . In particular we show, that a perfect squeezing over 
a finite solid angle leads to the same dynamics of the 
open system as an imperfect squeezing over the complete 
solid angle Air. (iii) A thermal mixture of coherent states. 
This example illustrates, that the equation of motion for 
the reduced probability distribution does not depend on 
the actual representation of the state of the environment. 



APPENDIX A: CALCULATION OF THE 
ENSEMBLE AVERAGE OF THE SECOND 
ORDER PROPAGATION OPERATOR 



2n to 



(A7) 



Thus, the constant N we have just defined is the mean 
number of photons in the modes with frequency uj s , where 
the mean is taken over the states of the ensemble and over 
the relative strength of the coupling. The imaginary part 
Si leads to the Stark shift. 

For the calculation of Gn we can use the commutator 
relations = Skk' and find: 



Cr U = C!,, - r£- I duj k h{uj k ) I dse^'-"^, (A8) 



where 



h(u)k) 



2ttuj 2 V 
(27rc) 3 7 



9k- 



(A9) 



Let us calculate first G22 '■ Combining Eq. ([46|) and ( |20| ) 
and performing the continuum limit yields 

G22 = -r^ J duj k n(u; k ) J ds e -*(— (Al) 

where we defined the mean number of photons n(to k ) with 
frequency to k as 



Thus we find 



n{u k ) 



2ttuj 2 V 
(27rc) 3 7 



[ cM gl{{blb k )) Ph[ . M , (A2) 



where k = (to k , k,X k ) and 

(sin 8 cos tp 
sin 9 sin tp 
cos 6 

and the mean quadratic coupling at frequency to s 
2<kuj 2 V 



(A3) 



7 



(27r C )3 ^ I g {u„k\\ h y 



(A4) 



Note that the above summations and integrations are 
restricted to modes k £ JCb- Since g\ ~ V^ 1 the expres- 
sions for n(to k ) and 7 are independent of the volume of 
quantization V. Using the identity 



Gu = -T(j(N + i)-i(S + S 1 )) (A10) 
where we used the definition of 7 (Eq. (|A4|)) and defined 

So =^p r d, k ^±^. (aid 

2tt to k 

The new constant So which appears here is the Lamb 
shift. 

In a similar way we can also calculate G12 = G\\ by 
defining 



m(u>) 



2ttuj 2 V 
(2ttc) 3 7 



]T / ' dn 9 l((b k bk)) Pk[ ., 0h (A12) 



where the summation and integration extend over all 
modes k £ /Cb, and 

M = i(, s ) + IP r du m{uJs+UJ \ (A13) 

7T ./ „ tO 



which leads to the result 



G12 = r-Me 2iw »*° 
G21 = r'~M*e~ 



(A14) 
(A15) 



2tt 



dukf{uk) / dre^ T 



5/(0) + ^ P / ''-v. /uv)/^. (A:,) 



where P denotes the Cauchy principal value we find 



G 2 2=-r(^N + iSi), 



(A6) 



with N — n(uj s ) and 
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FIG. 1. Correlation function {(a + o + (t)<j~ (t)ct~)) s for a 
coherently driven two level atom on resonance with Rabi fre- 
quency Q — IO7. The numerical result for 10 5 realizations 



of the simulation algorithm described in Sec. VA (errorbars) 
and analytical solution (thick line). 
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FIG. 2. Calculation of the first order correlation function 

{(cr + (r)(T _ )) s for a coherently driven two level atom on res- 
onance. This figure shows the CPU time in seconds vs. the 
relative error for the simulation algorith ms pro posed in Jl]] 
(Method I) and for our algorithm (cf. Sec. |VB2| , Method II). 
The solid lines represent the mean square deviation of the 
numerical solution from the exact solution and the dashed 
lines show the estimated standard deviation of the numerical 
solution. 
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FIG. 3. Resonance fluorescence spectra for a coherently driven two level atom with Rabi frequency Q = IO7 coupled to a 
squeezed vacuum for different squeezing parameters (mean photon number N, relative phase of squeezing and driving field 
4> = 2(<j> 3 — 4>l), and fraction e of solid angle which is squeezed): numerical results (solid line), and the vacuum spe ctrum (thin 
line). The numerical solution was obtained by simulating a stochastic process in a doubled Hilbert space (cf. Sec. V B jj ). 
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